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INTRODUCTION 

Adhesive  bonding  is  increasingly  used  to  fasten  structural 
components  together.  This  is  because  in  many  present  day  applications, 
conventional  fasteners  such  as  bolts,  rivets,  welds  etc.,  are 
unsuitable,  specially  if  the  components  are  made  of  polymeric  or 
composite  material.  Penetration  methods  (i.e.,  drilling  holes,  etc.) 
cause  high  stress  concentrations  and,  in  the  case  of  composites,  sever 
the  fiber  reinforcement  which  in  turn  reduces  the  strength  of  the 
joint.  On  the  other  hand,  bonded  joints  tend  to  be  damage-tolerant  due 
to  the  high  damping  behavior  of  the  adhesive  layer  and  less  expensive 
due  to  lower  fabrication  cost. 

Most  polymeric  adhesives  are  rate  sensitive  material  and  hence 
exhibit  viscoelasticity.  Furthermore,  certain  types  of  epoxy  resins 
have  been  found  to  be  nonlinearly  viscoelastic  in  character.  The 
nonlinear  viscoelastic  behavior  is  typified  by  a  stress-enhanced 
creep.  Basically,  at  elevated  stresses  the  material  moduli  seem  to 
soften  and  the  creep  progresses  at  accelerated  rates. 

At  present  there  are  numerous  computer  programs  available  for 
analyzing  bonded  joints.  However,  most  of  these  computer  codes 
incorporate  linearly  elastic  material  behavior,  and  some  allow  for 
nonlinearly  elastic  and  plastic  behavior.  Computer  programs  which 
incorporate  viscoelastic  material  behavior  are  quite  often  limited  to 
the  simple  spring-dashpot  type  of  model  for  linear  materials.  Such 


inaccurate  modelling  of  the  constitutive  behavior  of  the  structure  can 
seriously  compromise  the  accuracy  of  the  analytical  predictions. 

Ir  ‘■his  paper  we  are  conerned  with  the  development  of  a 
computational  procedure  which  is  capable  of  accurately  modelling  the 
linear  and  nonlinearly  viscoelastic  behavior  of  an  adhesive  layer  within 
a  bonded  joint.  The  adherends  can  be  modeled  as  a  linearly  elastic 
material  which  can  undergo  large  rotations  and  displacements.  The 
theory  and  computer  implementation  of  the  geometrically  nonlinear 
elastic  response  has  been  presented  in  a  previous  work  [1]  and  shall  not 
be  included  here.  The  rest  of  the  presentation  will  primarily  involve 
the  description  and  the  implementation  of  a  general  constitutive  law  for 
nonlinear  viscoelastic  materials. 

A  general  single  integral  constitutive  law  for  nonlinear 
viscoelastic  materials  systems  was  proposed  by  Schapery  [2].  The  law 
can  be  derived  from  fundamental  principles  using  the  concepts  of 
irreversible  thermodynamics.  A  comprehensive  review  of  the 
thermodynamics  basis  of  Schapery's  theory  has  been  presented  by  Hiel  et 
al.  [3]. 

The  present  study  deals  with  the  development  of  a  finite  element 
model  that  is  based  on  Schapery's  single  integral  constitutive  law. 
First,  a  stress  operator  that  defines  uniaxial  strain  as  a  function  of 
current  and  past  stress  is  developed.  Extension  to  multi-axial  stress 
state  is  accomplished  by  incorporating  Poisson's  effects,  resulting  in  a 
constitutive  matrix  that  consists  of  instantaneous  compliance,  Poisson's 
ratio  and  a  vector  of  hereditary  strains.  The  constitutive  equation 
thus  obtained  are  suitable  for  non-linear  finite  element  analysis. 

Plane  stress,  plane  strain  and  axisymmetric  formulations  are  included. 


A  brief  review  of  literature  is  presented  in  the  next  section  to  give 
the  background  for  the  present  study. 

BACKGROUND  LITERATURE 

An  analysis  of  adhesive  stresses  in  bonded  joints  was  first 
performed  by  Goland  and  Reissner  [4]  for  two  limiting  cases:  (i)  the 
case  in  which  the  adhesive  layer  is  so  thin  and  stiff  that  its 
deformation  can  be  neglected;  and  (ii)  the  case  in  which  the  adhesive 
layer  is  soft  and  flexible  and  the  joint  flexibility  is  mainly  due  to 
the  deformation  of  the  adhesive  layer.  In  the  first  case,  the  peel 
stress  is  found  to  be  very  high  at  the  edge  of  the  joint,  while  the 
shear  stress  is  zero.  In  the  second  case,  the  maximum  values  of  the 
peel  and  shear  stresses  occur  at  the  edges  of  the  joint.  The  Goland- 
Reissner  analysis  is  limited  to  identical  adherends,  the  joint-edge 
loads  are  not  in  equilibrium,  and  the  stresses  across  the  adhesive  layer 
are  constant. 

Erdogan  and  Ratwani  [5]  developed  a  one-dimensional  analytical 
model  for  calculating  stresses  in  a  stepped  lap  joint.  One  adherend  was 
assumed  to  consist  of  an  isotropic  material  and  the  second  of  an 
orthotropic  material.  Linear  elastic  conditions  for  the  materials  were 
assumed.  The  thickness  variation  of  the  stresses  in  both  the  adherends 
and  in  the  adhesive  was  neglected. 

Wooley  and  Carver  [6]  investigated  the  stress  distributions  in  a 
simple  lap  joint  using  the  finite  element  method.  They  assumed  that  the 
total  length  of  the  adherends  beyond  the  overlap  is  long  and  a  plane 
stress  condition  exists.  The  constant  strain  quadrilateral  element 
(obtained  by  combining  four  constant  strain  triangular  elements)  was 
used.  One  end  of  the  adherend  was  assumed  to  be  hinged  and  other  end 


was  allowed  to  move  freely  in  the  direction  parallel  to  the  original 
bond  line.  The  study  dealt  with  the  influence  of  Young's  moduli  ratios 
and  geometries  on  the  peel  and  shear  stress  distributions.  The  results 
compared  favorably  with  the  results  of  Goland  and  Reissner. 

Hart-Smith  (7)  improved  upon  the  approach  of  Goland  and  Reissner  by 
considering  a  third  free-body-diagram  for  the  adherend  outside  the  joint 
in  addition  to  the  two  free-body-diagrams  from  each  of  the  upper  and 
lower  halves  of  the  joint.  With  three  separate  sections  to  consider, 
three  relations  between  displacements  and  bending  moments  were 
obtained.  Additional  boundary  conditions  involving  displacements  and 
their  first  derivatives,  not  considered  by  Goland  and  Reissner,  were 
imposed  in  order  to  solve  for  the  additional  unknowns.  In  addition  to 
the  improvement  on  the  analysis  of  Goland  and  Reissner,  Hart-Smith  [7] 
also  estabished  the  quantitative  influence  of  adhesive  plasticity  in 
shear.  The  elastic-plastic  theory  used  by  Hart-Smith  predicts  an 
increase  in  joint  strength  and  was  shown  to  be  capable  of  explaining 
premature  failure  predictions  found  when  using  linear  elastic 
analyses.  The  quantitative  effects  of  stiffness  imbalance  were  also 
accounted  for. 

A  finite-element  stress  analysis  for  adhesive  lap  joints  using 
linear  elasticity  and  elasto-plasticity  theories  was  reported  by  Liu 
[8).  Stress  distributions  and  concentrations  in  the  adhesive  layer  for 
different  joint  parameters  (geometry,  material  properties,  and  loading 
conditions)  were  studied  and  compared. 

The  existence  of  stress  gradients  through  the  thickness  of  the 
adhesive  layer,  close  to  the  joint  edges,  was  observed  by  Adams  and 
Peppiatt  [ 9 ] .  They  subsequently  performed  a  linear  elastic  finite 


element  analysis  on  adhesively  bonded  lap  joints  employing  more  than  one 
element  through  the  thickness  of  the  adhesive  layer,  close  to  the  joint 
edges.  Adams  and  Peppiatt  [ 10 1  also  studied  the  adhesive  yielding  using 
an  iterative  elastic-plastic  finite  element  program  and  the  double  lap, 
bevel  and  scarf  joints.  The  adhesive  was  assumed  to  be  elastic- 
perfectly  plastic. 

A  nonlinear  analysis  of  single  and  double  lap  joints  was  presented 
by  Humphreys  and  Herakovich  [11]  using  the  finite  element  method.  The 
nonlinear  stress-strain  response  was  represented  by  a  Ramberg-Osgood 
approximation.  Mechanical  and  thermal  loadings  were  considered  but  only 
one  element  through  the  thickness  of  the  adhesive  layer  was  used. 

Allman  [12]  presented  an  elastic  stress  analysis  based  on  the 
strain  energy  density  of  a  particular  joint.  The  effects  of  bending, 
stretching  and  shearing  of  the  adherends  were  included,  and  the  shearing 
and  tearing  action  in  the  adhesive  was  accounted  for.  All  conditions  of 
stress  equilibrium  in  the  joint  and  stress-free  surface  conditions  were 
satisfied.  It  was  assumed,  however,  that  the  axial  stress  varies 
linearly  through  the  adherend  thicknesses  and  that  the  shear  stress  is 
constant  through  the  adhesive  thickness.  Allman  obtained  solutions  for 
the  single  lap  joint,  although  the  method  also  appears  to  be  applicable 
to  other  joint  configurations.  He  found  that  the  average  shear  stress 
concentration  is  11%  higher  than  that  of  Goland  and  Reissner's  first 
analysis,  while  the  average  peel  stress  at  the  joint  edge  is  67% 
lower.  Compared  with  the  second  analysis  of  Goland  and  Reissner, 
Allman's  method  yielded  a  shear  stress  concentration  of  15%  and  31%  less 
for  metal  and  composite  adherends,  respectively,  while  the  average  peel 
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stress  at  the  joint  edge  was  11%  higher  and  36%  lower  for  the  same  types 
of  adherends,  respectively. 

Phenomenological  considerations  were  discussed  by  Hart-Smith  [13] 
which  greatly  improve  our  understanding  of  the  sources  of  non-uniform 
load  transfer,  viz.,  adherend  extensitivity,  stiffness  imbalance  and 
thermal  mismatch.  He  also  explained  how  the  lightly  loaded  central  area 
of  the  joint,  away  from  the  joint  edges,  restricts  cumulative  creep 
damage,  and  suggests  that  this  region  is  vital  for  long  term 
durability.  The  amount  of  lightly  loaded  central  area  is  a  function  of 
the  overlap  length. 

Yuceoglu  and  Updike  [14]  developed  a  numerical  method  for 
determining  peel  and  shear  stresses  in  the  adhesive  of  double  lap, 
double  strap  and  stiffner  plate  joints.  Bending  and  transverse  shear 
were  included  in  the  analytical  model.  Shear  stresses  were  not  required 
to  drop  to  zero  at  the  joint  edges  after  reaching  peak  values  close  to 
the  edges.  Yuceoglu  and  Updike  maintained  that  an  analytical  model 
which  would  allow  the  shear  stresses  to  drop  to  zero  at  the  joint  edges 
would  give  approximately  the  same  or  slightly  lower  peak  values  of  shear 
and  peel  stresses.  Their  method  also  reveals  that  adherend  bending  has 
a  significant  effect  on  both  adhesive  shear  and  peel  stresses, 
especially  the  latter. 

Delale  and  Erdogan  [15,16]  performed  a  plate  analysis  simiar  to 
that  of  Goland  and  Reissner  on  the  single  lap  joint  assuming  linear 
elastic  adherends  and  a  linear  viscoelastic  adhesive.  Separate  stress 
distributions  were  calculated  for  membrane  loading,  bending,  and 
transverse  shear  loading.  They  further  extended  their  viscoelastic 


analysis  of  the  single  lap  joint  to  include  time-dependent  temperature 
variations. 

Gali  and  Ishai  [17)  performed  a  finite  element  analysis  on  a 
symmetric  doubler  model  with  linear  elastic  adherends  and  the  adhesive 
obeying  a  nonlinear  effective-stress-strain  relationship.  The 
effective-stress-strain  relationship  was  derived  from  stress-strain 
curves  obtained  by  tensile  and  shear  test  data,  and  based  on  the  Von 
Mises  deviatoric  energy  yield  criterion.  An  iteration  procedure  was 
applied  to  the  linearly  elastic  finite  element  problem  using  a  specific 
secant  modulus  for  each  element  separately.  The  secant  modulus  was 
found  from  the  corresponding  effective  strain  of  the  previous  solution 
and  the  corresponding  effective  stress  was  found  from  the  experimental 
stress-strain  curves.  Gal i  and  Ishai  analyzed  the  symmetric  doubler 
model  using  both  plane  stress  and  plane  strain  and  found  that  the  plane 
strain  solutions  converged  faster  and  yielded  less  conservative  results 
(i.e.,  lower  stresses)  than  the  plane  stress  solutions.  Nonlinear 
solutions  were  also  found  to  be  considerably  lower  than  the  linear 
solutions,  the  difference  being  more  pronounced  in  the  plane  stress 
case.  The  problem  was  also  solved  with  the  adhesive  following  an 
elastic-perfectly-plastic  effective-stress-strain  law.  The  difference 
between  these  results  and  those  of  the  continuous  nonlinear  effective- 
stress-strain  case  was  found  to  be  very  small. 

Nagaraja  and  Alwar  [18]  analyzed  a  tubular  lap  joint  with  the 
finite  element  method  assuming  linear  elastic  adherends  and  a  nonlinear 
biaxial  stress-strain  law  in  the  adhesive.  They  demonstrated  that  for 
low  stress  levels,  of  the  order  of  12%  of  the  fracture  stress,  the 
nonlinear  stresses  were  as  much  as  15%  lower  in  shear  and  8%  lower  in 


peel  than  the  linear  stresses.  Nagaraja  and  Alwar  [191  also  performed  a 
finite  element  analysis  on  a  single  lap  joint,  treating  the  adherends  as 
linear  elastic  materials  but  the  adhesive  as  a  linear  viscoelastic 
material.  The  relaxation  modulus  was  assumed  to  be  equal  to  the  inverse 
of  the  creep  compliance,  the  latter  being  obtained  experimentally. 

Francis  et  al.  [201  discussed  the  effects  of  a  viscoelastic 
adhesive  layer,  geometry,  mixed  mode  fracture  response,  mechanical  load 
history,  environmental  history  and  processing  variations  on  the  fracture 
processes  of  adhesively  bonded  joints.  However,  their  finite  element 
analysis  is  based  on  linear  elastic  fracture  mechanics. 

Dattaguru,  et  al.  [21]  studied  the  crack  lap  specimen  and  performed 
analyses  with  a  finite  element  program,  GAMNAS,  developed  in-house  at 
NASA-Langley.  The  program  GAMNAS  accounts  for  geometric  and  material 
nonl inearities  but  does  not  include  viscoelastic  capability.  Also, 
fracture  is  modeled  using  the  linear  elastic  fracture  mechanics  but  no 
failure  law  is  included. 

Botha,  Jones  and  Brinson  [22j,  Henriksen  [23],  Becker,  et  al.  [24], 
and  Yadagiri  and  Papi  Reddy  [25]  reported  results  of  the  viscoelastic, 
finite-element  analysis  of  adhesive  joints.  Henriksen  used  Schapery's 
[2]  nonlinear  viscoelastic  model  to  verify  the  experimental  results  of 
Peretz  and  Weitsman  [26]  for  an  adhesive  layer.  Becker  et  al.  [24] 
developed  a  finite-element  stress  analysis  program,  called  VISTA,  for 
adhesively  bonded  joints.  The  program  uses  the  nonlinear  viscoelastic 
model  of  Knauss  and  Emri  [27].  No  illustrative  or  validation  problems 
are  presented  to  demonstrate  the  two-dimensional  nonlinear  viscoelastic 
capabilities  of  VISTA.  The  works  of  Botha,  et  al.,  and  Yadagiri  and 
Papi  Reddy  are  limited  to  linear  viscoelastic  analysis. 


Pickett  and  Hollaway  [28|  presented  both  classical  and  finite 
element  solutions  for  elastic-plastic  adhesive  stress  distributions  in 
bonded  lap  joints.  Single,  double  and  tubular  lap  configurations  having 
both  similar  and  dissimilar  adherends  were  considered.  The  results  show 
how  the  development  of  adhesive  yielding  as  the  joints  are  loaded  to  a 
failure  condition.  The  detrimental  effect  of  adherend-st i f fness- 
imbalance  on  the  adhesive  stress  distributions  was  also  shown. 

An  approximate  method  to  analyze  viscoelastic  problems  has  been 
outlined  by  Schapery  [ 29 | .  In  this  method,  the  solution  to  a 
viscoelastic  problem  is  approximated  by  a  corresonding  elasticity 
solution  wherein  the  elastic  constants  have  been  replaced  by  time 
dependent  creep  or  relaxation  functions.  The  method  may  be  applied  to 
linear  as  well  as  nonlinear  problems.  Weitsman  [30]  used  Schapery's 
quasi-elastic  approximation  to  investigate  the  effects  of  nonlinear 
viscoelasticity  on  load  transfer  in  a  symmetric  double  lap  joint. 
Introducing  a  stress-dependent  shift  factor,  he  observed  that  the 
enhanced  creep  causes  shear  stress  relief  near  the  edges  of  the  adhesive 
joint. 

Schaffer  and  Adams  [311  carried  out  a  nonlinear  viscoelastic 
analysis  of  a  unidirectional  composite  laminate  using  the  finite  element 
method.  The  nonlinear  viscoelastic  constitutive  law  proposed  by 
Schapery  [ 2 ]  was  used  in  conjunction  with  el astopl ast i c  constitutive 
relations  to  model  the  composite  response  beyond  the  elastic  limit. 

Ghoneim  and  Yu  Chen  [ 32 |  developed  a  viscoelastic-viscoplastic  law 
based  on  the  assumption  that  the  total  strain  rate  tensor  can  be 
decomposed  into  a  viscoelastic  and  a  viscoplastic  component.  A  linear 
viscoelasticity  model  is  used  in  conjunction  with  a  modified  plasticity 


model  in  which  hardening  is  assumed  to  be  a  function  of  viscoplastic 
strains  as  well  as  the  total  strain  rate.  The  resulting  finite  element 
algorithm  is  then  used  to  analyze  the  strain  rate  and  pressure  effects 
on  the  mechanical  behavior  of  a  viscoelastic-viscoplastic  material. 

Analysis  of  crack  growth  in  viscoelastic  media  are  mainly  limited 
to  linear  isotropic,  homogeneous  materials.  Schapery  [331  proposed  the 
use  of  parameters  similar  to  the  J  integral  for  quasi-static  crack 
growth  in  a  class  of  nonlinear  viscoelastic  materials  subject  to  finite 
strains. 

Czarnocki  and  Piekarski  ( 34 1  used  a  nonlinear  elastic  stress-strain 
law  for  three-dimensional  failure  analysis  of  a  symmetric  lap  joint. 
Taking  into  account  the  variation  of  Poisson's  ratio  with  strain  within 
the  adhesive,  the  authors  concluded  that  the  failure  of  the  adhesive 
layer  originates  in  the  central  plane  of  a  joint  (at  the  front  edge). 

It  was  also  observed  that  the  joint  width  does  not  have  any  effect  on 
the  stress  peaks  in  the  central  plane  and  that  the  application  of  a 
weaker  but  more  flexible  adhesive  results  in  higher  load  carrying 
capacity  and  lower  stress  concentrations  in  the  adherends. 

NONLINEAR  VISCOELASTIC  FORMULATION 

The  uniaxial  single  integral  constitutive  equation  can  be  stated 
as, 

t 

c(t)  =  g0(t)o(t)Do  +  gi(t)  /  Dc(/  -  US)  5-  [g2(s)o(s)]dS  (1) 

In  Eq.  (1),  e(t)  represents  uniaxial  kinematic  strain  at  current  time 
t,  o(t)  is  the  Cauchy  stress  at  time  t,  0Q  is  the  elastic  compliance  and 
0(h))  is  a  transient  creep  compliance  function.  The  factor  g  (t) 
defines  stress  and  temperature  effects  on  elastic  compliance  and  is  a 
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measure  of  state  dependent  reduction  (or  increase)  in  stiffness. 


Transient  (or  creep)  compliance  factor  g^(t)  has  similar  meaning, 
operating  on  the  creep  compliance  component.  The  factor  g^(s)  accounts 
for  the  influence  of  load  rate  on  creep,  and  depends  on  stress  and 
temperature.  The  functions  ^  and  wS  are  reduced  times  at  t  and  s, 
respectively,  and  are  defined  by, 

^  =  Jrt  (a*  )_1ds  (2) 

0  °T 

where,  asT  is  a  time  'shift  factor'.  This  function  modifies 
viscoelastic  response  as  a  function  of  temperature  and  stress. 
Mathematically,  a^  shifts  the  creep  data  parallel  to  the  time  axis 
relative  to  a  master  curve  for  creep  strain  versus  time. 

The  constitutive  law  in  Eq.  (1)  can  be  expressed  in  the  following 
operational  form: 

E  =  F(o)  (3) 

where  F(a)  is  a  stress  operator,  defined  by, 

F(o)  =  Ojo  +  E  (4) 

The  detailed  derivation  of  Eq.  (4)  from  Eq.  (1)  is  given  below. 

The  transient  creep  compliance,  Dc(v),  can  be  expressed  in  the 
following  exponential  form, 

Dr(«0  =  l  D  11  -  e  X?>]  (5) 

r 

Substituting  (5)  in  (1),  gives, 

t  (u^-|iS)  . 

e  =  90Do°  +  gl  /  lDr[l~e  I  [g2(s)o(s) Ids  (6) 

0  r 

Letting  the  product  g^a  be  expressed  as  G  and  simplifying  the 
integrand  on  the  right  hand  side  of  Eq.  (6)  yields, 
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t.  t  -V  ( ^  -*))  -,  . 

'  -  30V  ♦  *1  I  °r  -0  3s  G<S>dS  •  31  l  °r  e  1f5i  d5 


The  third  integration  term  on  the  right  hand  side  of  Eq.  (7)  is 
now  separated  into  two  parts,  the  first  part  having  limits  from  0  to 
to  (t  -  At)  and  the  second  integral  spanning  only  the  current  load  step. 


i.e.,  from  (t  -  At)  to  t.  Hence, 

*  ~Vr(i,t  -  dG(s) 
J  e  -hH 


_t-At  -Vr(/-,S)  dQ(s) 


-xr(»  ~v  )  dGfs) 

♦  .  e  ds 
t-At  05 


Now,  the  first  term  on  the  right  hand  side  of  Eq.  (8)  can  be  rewritten 


*1*1  ds 


f  .  t  s 

-t-'t  "V  V  dG(s) 
e  e  1  ds 

n  dS 


f  ,  t  ,  t-At  (  S  ,  t-At 

-t_at  ‘V’  V  V"  'V  dG(s)  • 
•e  e  e  e  ds 

0  ds 


,  t  t-Ats  .  ,  t-At  S> 
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t-At  =  At"^S)  dGisl 

H  y  y 

r  0 


ds 


ds 


(ID 


The  second  integral  on  the  right  hand  side  of  Eq.  (8)  is  now 
integrated  by  parts.  To  carry  out  the  integration,  it  is  assumed 
that  G  varies  linearly  over  the  current  time  step  At.  Hence, 


I  e 

t-At 


dGLslds 

ds 


-A  (/-^S) 

.  dG(s)  e  _ 

ds  A_ 


t  t  2  ~xr 

_  J  dZG(s)  e 

t-At  t-At  ds2  Vr 


ds 


=  dG(t)  1_  _  dG(t-At)  e 


x  ,t-At. 

-Ar(4>  -il>  ) 


dt  a  dt 


-A  A4) 

.  dGiti  |1  -  e  | 


dt 


(12) 


In  arriving  at  the  second  step,  the  fact  that  G(s)  is  assumed  to  be 
linear,  hence  its  second  derivative  is  zero,  is  used.  Since  G(t)  has 
been  assumed  to  be  a  linear  function  of  time  over  the  current  load- 
step  At,  and  the  reduced  time  v  is  proportional  to  t,  we  can  write, 

dG(t)  _  G(t)  -  G(t-At) 

dt  ,t  t-At 

il>  - 


or. 


dG(t)  =  G(t)  -  G(t-At) 


dt 


(13) 


Ail) 
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Substituting  (13)  into  (12),  gives. 


t  -x  UV) 


■XrAil) 


f  e  r'T  T  §  ds  =  [ G ( t )  -  G(t-At)][^---  e- - ] 


t-At 


X  A4; 

r 


or. 


t  -x 
J  e 


§  ds  =  [ G ( t )  -  G(t-At) 18^  (14) 


t-At 


ds 


where 


6 


t  _  1  -  e 


-X  Ait; 


r  XrA^ 


(15) 


Substituting  (11)  and  (14)  back  into  (7)  and  writing  G  =  g  0 


e  =  goDo°  +  gl  l  Drg2° 

-X  Ail)^  . 

-  9!  I  Dr{e  r  q^At  +  [g2(t)o(t)  -  g2(t-At)a(t-At) ]ej}  (16) 

Collecting  those  terms  in  Eq.  (16)  that  are  multiplied  by  current 
stress  a  yields. 


e  =  [go°o  +  glg2  I  Dr  -  glg2  l  DrBr ]o 

r  r  t 

+  91^  Dr(92^t_At)Bra(t'at)  ~  e  qr _AtH  (17) 

Defining  instantaneous  compliance  Dj  as  the  compliance  term 
multiplying  the  instantaneous  stress  o,  and  the  remaining  terms  in  Eq. 
(17)  as  hereditary  strains  E,  yields. 


e  =  Djo  +  E  (18) 

where, 

°I  ■  9o°o  +  9192  l  °r  -  3^2  l  OrsJ 
r  r 


(19) 


if  ^  ^  ^  v.  *  :.v.  v  \^v,v.^viv.^viv,VAiyA^v.\v.t.v 


i 


4.  -X  A*  f  ... 

E  =  g^I  Dr(g2(t-At)s;o(t-At)  -  e  qj"4t]} 

r 


(20) 


Hence,  Eq.  (18)  expresses  Schapery's  single  integral  constitutive 
law  in  terms  of  a  stress  operator  that  includes  instantaneous  compliance 
and  hereditary  strains. 

.th 


.t-At  • 


It  is  to  be  noted  that  the  term  q^;”  in  Eq.  (20)  is  the  rT 
component  of  the  hereditary  integral  series  at  the  end  of  the  previous 
load  step  (i.e.  at  time  equals  t  -  At).  The  expression  for  the 
hereditary  integral  at  the  end  of  the  current  load  step  (i.e.  at  time  t) 
can  be  derived  in  the  form  of  a  recurrence  formula  as  shown  below. 


By  definition  [see.  Eq.  (11)], 


.  t  -x  (i|)  -tj)S)  ,r 

<L  ■  J  e  —  ds 

r  0 


ds 


t-At  -x>V)  p  t  -x  (^-id5)  Hr 

_  r  „  r  dG  ,  ,  r  _  r'  dG  . 

-  J  e  ds  ds  J  e  ds  ds 

0  as  t-At  ds 


Using  the  results  from  Eqs.  (11)  and  (14),  the  above  equation  reduces 
to. 


“X  Ail) 


qj:  =  e  r  qJ'At  +  [g2(t)o1j(t)  -  g2(t-At)o(t-At)leJ;  (21) 


t  ... 


where  Br  is  defined  by  Eq.  (15). 


MULTIAXIAL  STRESS  FORMULATION 


As  mentioned  earlier,  the  constitutive  law  derived  in  the  preceding 
section  holds  true  for  uniaxial  state  of  stress.  In  order  to  formulate 
a  stress  strain  relationship  for  a  multiaxial  stress  state  each  strain 
component  is  assumed  to  be  a  linear  function  of  the  stress  operators. 
Therefore,  as  in  linear  elastic  analysis,  Poisson's  effect  is 
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incorporated.  Hence,  the  multiaxial  stress  strain  law  is  fully  defined 
by  the  matrix  relationship, 

{e}  =  Dj[Hl{o}  +  [N] {E}  (22) 

In  Eq.  (22),  {e}  is  a  vector  containing  the  algebraic  difference  of 
kinematic  strains  {e}  and  thermal  strains  { & . ^ e } , 

® )  «  (e22  ®)  *  •  (e33  -  9)}  (23) 

while,  {o}  contains  four  components  of  Cauchy  stress, 

{a}  =  {°u»  °22’  J12*  °33 ^  (24) 

and  {E}  is  a  vector  of  hereditary  strains,  components  of  which  are 
defined  by  the  equation, 

(E}T  =  [Eu,  E ^  12 *  E33)  <25) 

Note  that  all  quantities  are  functions  of  current  time,  t. 

The  matrix  [N]  is  a  4x4  matrix  given  by. 


■  1  -V  0  -M  • 

[N]  =  1  0  _v  (26) 

0  0  2( I+v)  0 

■  -v  -v  0  1  - 

where  v  =  v(t)  is  the  Poisson's  ratio  at  time  t. 

It  is  to  be  noted  that  the  definitions  (22)  through  (26) 
incorporate  possible  states  of  plane  stress,  plane  strain  and  rotational 
symmetry. 


ir^7vrir*Tr7*TrT*Trr7Tm  »?t  'rm>' 7  v"  •-  t.^.v  v  wr,  v  v  v  x-^v vti 


Inversion  of  Eq.  (22)  gives  a  constitutive  relationship  which, 
written  in  a  matrix  form,  is. 


M  =  (Ml {e}  -  {E} 

0T 


(27) 


where  [ M 1  is. 


(Ml  = 


Dj(1+m) (1-2v) 


1-v 

V 

V 

1- 

0 

0 

-  V 

V 

0 

0 


0 


0 

1-v 


(28) 


Equations  (27)  and  (28)  provide  a  general  constitutive  law  and  can 
be  applied  to  either  plane  stress,  plane  strain  or  axisymmetric 
problems.  For  plane  strain,  kinematic  strain  component  e  —  is 
identically  zero.  Hence,  corresponding  stresses  may  be  computed  by 
setting  =  0.  Since  for  plane  stress,  o33  is  identically  zero,  the 
kinematic  strain  can  then  be  evaluated  from  Eq.  (27)  and  (28)  as, 

=  (1  +  v)(l  -  2v)  r  _  v  /  .  e)  -  "  (r  9)  +  0 

b33  1  -  v  '  11  e)  l-o  (e22  9 

(29) 
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(1  -  vL) 


FINITE  ELEMENT  SOLUTION  ALGORITHM 


Basic  Formulation 


The  principle  of  virtual  work  states  that  for  a  system  to  be  in 
equilibrium  (see  ( 1,351) , 


6W.  .  =  6W  .  . 

internal  external 


or. 


fiW.  .  -  «W  =  0 
int  ext 


(30) 
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* 


V 
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For  any  elastic  structure, 

5Wint  =  $  o6edv 


(31) 


where,  6e  is  the  variation  in  the  strain  caused  by  the  virtual 
displacement  au. 

Equation  (31)  can  be  written  in  matrix  form  as. 


6U 


int  =  -T 


(32) 


The  strain-displacement  relations  can  be  expressed  as, 

ifie}  =  [B|(6u[  (33) 


where  [B]  is  the  transformation  matrix  relating  strains  to  displacements 


r_—  0 

ax 


[B]  = 


a_  a_ 

ay  ax. 

Needless  to  mention,  the  above  relationship  holds  only  for 


0  f- 

ay 


(34) 


geometrically  linear  strain-displacement  equations.  Substituting  Eq. 
(33)  in  Eq.  (32)  gives, 

6Wint  =  { 5U }  T[  B  ] T  {o }  dv  (35) 

The  variation  in  external  work  due  to  virtual  displacement  { su }  is 
given  by. 


6W 


ext 


I  f«u] 


■Fext^dv 


(36) 


Substituting  appropriate  finite-element  interpolation  of  the 
displacements  (see  (361)  into  Eqs.  (35)  and  (36),  and  substituting  the 
result  into  Eq.  (30),  we  obtain  the  finite-element  equations. 


!«u|T(;  [BlVldv  -  {FextU  -  (0] 


(37) 


where  { 6u }  is  a  vector  of  virtual  displacements,  and  {frextl  Is  a  vector 
of  externally  applied  forces  at  the  boundary. 

If  the  stress  vector  {o}  in  Eq.  (37)  is  now  replaced  by  the  multi- 
axial  constitutive  law,  Eq.  (27)  derived  previously,  then  Eq.  (37)  can 
be  expanded  as, 

{ <5 u } T { _f  [ B ] T[M]  {e }dv  -  /  [B]T[M]  {9}dv 
v  v 

-  J  ±-  [ B ] T { E } d v  -  (F  }}  =  {0}  (38) 

v  D.  ext 

or  simply,  1 

mVi  -  (o) 

where 

{R}  =  J  l B 1 T { o}dv  -  {Fext}  (40) 

is  the  residual  force  vector. 

Inspection  of  Eq.  (36)  reveals  that  it  is  a  function  of  the 
instantaneous  compliance  Dj  at  the  current  time  step.  But  from  Eq.  (19) 
it  is  observed  that  Dj  is  a  function  of  gQ,  g^,  g^,  and  the  reduced  time 
and  therefore  dependent  on  the  stress  and  temperature  at  the  current 
time  step.  Hence,  Eq.  (38)  is  nonlinear  and  cannot  be  solved  by  the 
direct  method.  The  solution  must  be  obtained  by  means  of  an  iterative 
method  which  is  presented  in  the  following  section. 

Iterative  Scheme 

In  determining  the  solution  to  Eq.  (38)  at  any  time  step  t,  the 
following  approach  is  adopted.  If  {R.+^}  is  the  vector  of  unbalanced 
forces  at  iteration  i  +  1  for  time  t,  then, 

{Ri+1}  =  {Ri 1  +  [KT]{aui} 


19 


where  [Kyi  is  the  tangent  stiffness  matrix, 


(K  I  =  (42) 

d{u} 

in  which  {u}  is  a  vector  of  nodal  displacements  at  time  t.  It  can  be 
shown  from  Eq.  (38)  that  [Ky]  is  given  by, 

(Kyi  =  J  [B]T[M][B]dv  (43) 

v 

Since  we  seek  for  a  solution  in  which  {R^}  =  {0},  incremental 
displacements  { au ^ }  are  evaluated  from, 

{AUi}  =  -  [ Ky ] - 1 { R. }  (44) 

The  displacement  is  updated  after  each  iteration  by, 

{ui+i}  =  {u.}  +  {au. }  (45) 

As  can  be  perceived,  the  above  solution  scheme  is  essentially  the 
Newton-Raphson  iterative  solution  algorithm. 


Solution  Algorithm 


The  complete  solution  procedure  for  each  individual  time  step  is 
presented  in  a  logical  step-wise  fashion  and  can  be  used  directly  for 


programming  purposes: 

1.  At  the  beginning  of  each  time  step,  the  stress  vector 

{a}  from  the  previous  time  step  is  accessed.  Note  that  for 
the  initial  or  starting  time  step,  the  stress  vector 
a(t  -  At)  denotes  the  initial  stress  state  at  t  =  0,  given 
by  {o0}.  Since  it  is  customary  to  assume  a  stress  free  state 


to  exist  at  the  start  of  the  solution,  {o0}  is  usually  set  to 
be  zero. 

2.  Temperature  T  at  time  t  is  computed  from  T  =  f(t)  which  is 
supplied  by  the  user  for  problems  involving  thermal  loads. 

3.  The  parameters  gQ(t),  g ^ ( t ) ,  g^t),  and  a^,  which  are  known 
functions  of  temperature  and  stress,  are  evaluated  next,  using 
the  stress  vector  obtained  from  previous  load  step. 

4.  Assuming  a^  to  be  a  linear  function  of  time  over  the  time 

step  At,  the  average  value  of  shift  factor  is  given  as 

a^-r  =  (atZAt  +  atT)/2  and  the  change  in  reduced  time 
oTavg  °T  t  “T  t 

is  computed  as  =  At/a  T  .  In  order  for  this  assumption 

0 '  avg 

to  be  valid  At  should  be  made  sufficiently  small. 

5.  Hereditary  integral  {q£}  is  computed  using  the  recursive 
formula  given  by  Eq.  (21). 

6.  (F  .}  =  x [ F  }  where  x  is  the  load  factor  that  corresponds 

L  extJ  appJ 

to  the  time  step  under  consideration. 

7.  The  residual  vector  {R}  is  computed  for  each  element  as, 

We  ■  {Fext}e  -  J  I B 1 T { o}edv 

e 

8.  The  tangent  stiffness  matrix  [ ] e  =  J  [B]T[M] [B]dv. 

e 

9.  Incremental  displacement  { au }  =  [K^I”^{R}. 

10.  Total  displacement  {u } .  =  {u}.  ^  +  {au}..  where  the  subscript  i 
denotes  the  number  of  iterations. 

11.  The  strains  and  stresses  are  computed  using  the  known 
displacement. 

II  {au  }n 

12.  Steps  3  through  12  are  repeated  till  -  <  Tolerance. 

il{u.}n 


13.  Solution  proceeds  to  the  next  time  step  for  which  steps  1 
through  12  are  repeated. 

SAMPLE  PROBLEMS 

Validation  of  Linear  Viscoelastic  Model 

The  nonlinear  constitutive  law  due  to  Schapery  may  be  linearized  by 
assuming  that  the  nonlinearizing  parameters  gQ,  g^,  and  g 2  have  a  value 
of  unity.  In  addition,  the  stress  dependent  part  of  the  exponent  in  the 
definition  of  the  shift  function  is  set  to  zero.  Consequently,  the 
constitutive  law  reduces  to  the  hereditary  integral  form  commonly  used 
to  describe  a  linear  viscoelastic  material. 

Two  test  cases  are  used  to  validate  the  linear  viscoelastic 
analysis  capability  implemented  in  the  present  finite  element  program 
named  NOVA.  In  the  first  case,  the  tensile  creep  strain  in  a  single 
eight  noded  quadrilateral  element  was  computed  for  both  the  plane  stress 
and  plane  strain  cases  using  the  program  NOVA.  The  results  were  then 
compared  to  the  analytical  solution  for  the  plane  strain  case  presented 
in  [ 37 1 .  A  uniform  uniaxial  tensile  load  of  13.79  MPa  was  applied  on 
the  test  specimen.  A  three-parameter  solid  model  was  used  to  represent 
the  tensile  compliance  of  the  adhesive.  The  Poisson's  ratio  was  assumed 
to  remain  constant  with  time.  The  following  time  dependent  functions 
were  used  in  [ 37 ]  (see  [ 38 1  for  additional  experimental  data  on  FM- 73) 
to  represent  the  tensile  compliance  and  the  Poisson's  ratio  for  FM-73M 

nm  -  G°  ,  f  Gl  in  r-t/0.85, 

D(t)  '  2[l+v(t)]  +  Mi+v(t)]}(1  ‘  6  ) 


( 3K(t 

L-‘] 

2G(t 

,3K  t 

1  rrt\ 

L  +  11 

at  72°C: 


where  G(t)  and  K(t)  are  the  shear  and  bulk  modulus  (mm/mm/MPA) 
respectively.  The  analytical  solution  to  the  creep  problem  for  the 
plane  strain  case  is  given  in  [37]  as: 
e (t)  =  2.728  x  10'2  +  1.334  x  10"2  e"t/0>85  -  2.659  x  10"4  e~t/0‘3921 
It  is  to  be  noted  that  for  the  three-parameter  solid  charac¬ 
terization  of  FM-73M  the  value  of  the  Poisson's  ratio  increases 
significantly  with  :ime.  Consequently,  the  prony  series  coefficients 
for  the  tensile  compliance  also  change  with  time.  In  the  present 
formulation  compliance  coefficients  are  assumed  to  be  independent  of 
time.  Hence  two  discrete  values  of  the  Poisson's  ratio  are  used  to 
match  the  exact  solution  for  few  initial  time  steps  and  final  time 
steps.  The  values  of  the  Poisson's  ratio  chosen  for  this  purpose 

are  v  =  Lim  v(t)  =  0.32  and  v  =  Lim  v(t)  =  0.417.  Figure  la  shows 
0  t-0  *  t~ 

the  creep  curve  for  v  =  0.417  for  both  plane  strain  and  plane  stress 
finite-element  analyses.  As  expected,  the  plane  strain  results  exhibit 
close  agreement  with  the  exact  solution  for  large  values  of  time, 
followed  by  progressive  deteriorat ion  of  predicted  value  as  one  moves 
towards  smaller  values  of  time.  The  finite  element  results  for  the 
plane  stress  case  points  to  the  fact  that  the  strains  are  higher  for 
plane  stress  than  for  plane  strain. 

Figure  lb  shows  the  creep  curve  corresponding  to  v  =  0.32.  In  this 
case  the  finite  element  predictions  are  accurate  only  for  first  few  time 
steps  and  deviates  more  and  more  from  the  analytical  solution  as  time 
increases.  This  is  not  surprising  since  the  choice  of  Poisson's  ratio 
for  this  case  makes  the  comparison  meaningful  only  when  t  is  small. 

The  above  results  indicate  that  the  program  NOVA  provides 
reasonably  accurate  results  in  regions  where  the  input  parameters  are 


23 


accurate,  and  that  the  variation  of  Poisson's  ratio  during  the  period  of 
analysis  may  cause  significant  deviations  from  the  actual  solution. 

Next,  the  Model  Joint  analysis  problem  presented  in  [37)  was  used 
as  the  second  validation  example.  In  this  case,  a  linear  viscoelastic 
finite  element  analysis  was  carried  out  on  a  model  joint  under  a 
constant  applied  load  of  4448  N  giving  an  average  adhesive  shear  stress 
of  13.79  MPa.  The  specimen  geometry,  discretization  and  boundary 
conditions  are  shown  in  Fig.  2.  The  thickness  of  the  adhesive  layer  is 
taken  to  be  0.254  mm.  A  nine  parameter  solid  model  was  used  to 
represent  the  tensile  creep  compliance  of  FM-73  at  72  C  and  is  given  by: 

D(t)  =  0.5988  x  10"3  +  1.637  x  10"5  (1  -  e‘t/0-0L) 

+  0.6031  x  10'4  (1  -  e't/0,1) 

+  0.9108  x  10~4  (1  -  e*171,0) 

+  2.6177  x  10'4  (1  -  e't/10-°) 

The  Poisson's  ratio  is  assumed  to  have  a  value  of  0.417  and  remains 

constant  with  time. 

Figures  3  and  4  contain  plots  of  the  bond  normal  and  shear 
stresses,  respectively  for  t  =  50  secs,  and  t  =  60  min.  of  loading. 

These  stresses  represent  the  value  at  1/16  the  thickness  from  the  upper 
adhesive  adherend  interface.  The  sharp  peak  at  the  left  hand  edge  is 
due  to  the  singularity  caused  by  the  presence  of  a  re-entrant  corner  in 
the  vicinity  of  the  edge.  These  results  are  in  good  agreement  with  the 
results  presented  in  [37]  which  uses  the  linear  viscoelastic  finite 
element  eode,  MARC. 


Validation  of  Nonlinear  Viscoelastic  Model 

In  order  to  validate  the  nonlinear  viscoelastic  model,  three 
uniaxial  test  cases  are  analyzed.  The  results  are  compared  with  the 
laboratory  tests  conducted  on  similar  specimens  by  Peretz  and  Weitsman 
[26].  The  material  properties  used  in  the  verification  analysis  are 
those  reported  in  [23 [.  The  creep  data,  together  with  other  relevant 
material  properties,  are  given  in  Table  1.  A  constant  value  for  the 
Poisson  ratio  is  assumed  for  the  adhesive.  The  results  from  a  linear 
viscoelastic  analysis  are  also  presented  for  comparison. 

In  the  first  verification  test,  a  uniaxial  stress  of  10  MPa  is 
applied  to  the  adhesive  coupon  for  1200  secs.,  followed  by  a  step 
increase  to  26.6  MPa  for  a  further  1200  secs.  The  temperature  of  the 
specimen  is  held  constant  at  50°C  and  is  assumed  to  be  uniform 
everywhere.  The  finite  element  predictions  for  this  test  are  plotted 
together  with  the  experimental  data  in  Fig.  5.  The  predictions  are  in 
good  agreement  with  the  experimental  results  of  Peretz  and  Weitsman  [261 . 

The  second  test  involves  creep  predictions  under  simultaneously 
varying  stress  and  temperature,  both  increasing  linearly  with  time.  The 
temperature  is  again  assumed  to  be  uniform  throughout  the  test 
specimen.  The  finite  element  predictions  (linear  and  nonlinear)  and 
experimental  data  are  compared  in  Fig.  6.  There  is  a  good  agreement 
between  the  two  sets  of  data. 

The  third  test  involves  creep  under  a  constant  stress  of  10  MPa 
with  a  linearly  varying  temperature  as  a  function  of  time.  Figure  7 
shows  the  strain  vs.  time  curves  obtained  in  the  experiments  and  finite 
element  analysis.  Satisfactory  agreement  between  the  experimental 
results  and  the  analysis  is  observed. 
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After  the  successful  validation  of  linear  and  nonlinear 
viscoelastic  models,  the  program  is  used  to  investigate  stresses  and 
strains  in  a  bonded  cantilever  plate  undergoing  bending  due  to  a 
uniformly  distributed  transverse  load  (see  in  Fig.  8).  The  aluminum 
plates  are  bonded  together  by  means  of  a  thin  layer  of  adhesive.  The 
material  properties  used  for  aluminum  are  given  in  Table  2.  The  creep 
data  used  for  the  adhesive  is  the  same  as  the  one  employed  for  the 
verification  tests.  The  finite  element  mesh  used  for  the  problem  is 
also  shown  in  Fig.  8.  A  non-uniform  mesh  is  used  in  anticipation  of  a 
large  stress  gradient  near  the  clamped  end.  Eight-noded  quadri lateral 
isoparametric  plane  strain  elements  are  used  for  the  analysis.  The 
aluminum  plates  are  assumed  to  be  linearly  elastic  but  undergoing  large 
displacements,  and  the  adhesive  is  assumed  to  be  a  nonlinearly 
viscoelastic  material.  Both  the  applied  load  and  the  temperature  are 
assumed  to  be  constant  with  time. 

Figures  9  to  12  depict  the  results  of  this  analysis  over  a  period 
of  600  secs.  Figure  9  contains  plots  of  axial  stresses,  in  aluminum  and 
adhesive,  plotted  along  the  beam  axis.  Figures  10  and  11  contain  the 
variation  of  the  shear  stresses  and  strains,  respectively,  along  the 
axis  of  the  beam.  It  can  be  seen  that,  the  strains  in  the  adhesive  are 
an  order  of  magnitude  larger  than  the  strains  at  a  corresponding  axial 
location  in  the  aluminum.  Furthermore,  these  strains  increase  in 
magnitude  with  an  increase  in  time.  While  the  stresses  and  strains 
within  the  aluminum  plate  (approximately)  conform  to  the  pred't’ons  ;‘ 
elementary  beam  theory,  the  corresponding  strains  and  stresses 
the  adhesive  layer  exhibit  significantly  different  beha^'or.  'v 


stress  within  the  adhesive  change  from  compressive  to  tensile  along  the 
beam  axis  away  from  the  clamped  end.  In  the  neighborhood  of  the  free 
edge,  the  axial  stress  reaches  a  peak,  followed  by  a  sharp  drop  to  zero 
at  the  free  edge.  The  axial  and  shear  stresses  within  the  adhesive 
exhibit  a  certain  amount  of  relaxation  as  time  progresses. 

Figure  12  shows  through-the-thickness  variation  of  the  shear  stress 
at  two  different  locations,  x  =  13  mm  and  89  mm,  along  the  beam  axis. 

The  area  under  each  curve  is  obtained  by  means  of  numerical  integration 
to  verify  the  equilibrium.  The  total  vertical  force  at  each  section  did 
not  change  with  time  and  was  reasonably  close  to  the  actual  applied 
force  in  the  transverse  direction.  It  is  clear  from  the  results  that 
the  adhesive  layer  experiences  the  largest  shear  stress  close  to  the 
free  end  than  near  the  clamped  end. 

Linear  and  Nonlinear  Analysis  of  a  Model  Joint: 

The  loading,  boundary  conditions  and  specimen  geomtry  used  in  this 
analysis  is  the  same  as  the  one  used  in  the  earlier  model  joint  (see 
Fig.  1).  In  addition,  the  same  nine  parameter  solid  model  was  used  in 
this  analysis.  A  linear  viscoelastic  finite  element  analysis  was 
carried  out  over  a  period  of  one  hour  at  a  constant  applied  load  of 
3336  N.  The  results  for  the  linear  analysis  are  shown  in  Figs.  13-14. 
The  sharp  peak  at  the  left  hand  edge  is  due  to  the  singularity  caused  by 
the  presence  of  a  re-entrant  corner.  All  stress  plots  show  the  same 
basic  t^end  in  that  the  stresses  are  attempting  to  redistribute 
themselves  to  achieve  a  more  uniform  distribution. 


For  the  nonlinear  viscoelastic  analysis  of  the  model  joint,  the 
same  specimen  geometry  and  material  properties  were  employed.  However, 
the  non l inearizing  parameters  and  the  shift  function  were  no  longer  held 
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constant,  but  were  allowed  to  change  with  the  current  stress  state 
within  the  adhesive  layer.  The  results  from  this  analysis  are  presented 
in  Figs.  15-16.  It  is  immediately  apparent  that  the  effect  of  the 
nonlinearity  causes  a  'softening'  of  the  adhesive,  leading  to  a  response 
that  is  less  stiff  compared  to  the  linear  case.  Hence,  even  though  the 
applied  load  is  the  same,  the  shearing  strain  for  the  nonlinear  case  is 
significantly  larger  as  compared  to  the  linear  case  (Figs.  14  and  16). 
Moreover,  the  increment  in  creep  strain  for  the  nonlinear  case  is  0.0058 
as  compared  to  0.0041  for  the  linear  case  over  the  same  period  of 
time.  This  is  exactly  what  is  expected  since  the  nonlinear  model  takes 
into  account  the  acceleration  of  creep  caused  by  the  stresses  within  the 
adhesive. 

The  effect  of  the  nonlinearity  on  the  stress  curves  (Figs.  13  and 
15)  is  to  create  a  more  uniform  stress  distribution  by  reducing  the 
stress  peaks  near  the  edges  while  increasing  the  stresses  at  the  mid¬ 
section  of  the  overlap.  The  significant  reduction  of  the  stress  peaks 
effected  by  the  nonlinear  model  is  very  important  from  a  design  point  of 
view  since  the  reduction  of  stress  levels  at  the  critically  stressed 
regions  results  in  an  improved  joint  efficiency. 

Stress  Analysis  of  a  Thick  Adherend  Lap  Shear  Specimen 

The  geometry,  boundary  conditions  and  finite  element  mesh  of  the 
thick  adherend  lap  shear  specimen  are  shown  in  Fig.  17.  The  load  and 
the  temperature  are  kept  constant  with  time.  The  material  used  for  the 
adherend  is  aluminum,  and  the  adhesive  is  FM-73.  The  properties  of 
these  materials  have  been  1.  ted  in  Tables  1  and  2.  The  aluminum  is 
assumed  to  be  linearly  elastic  and  the  adhesive  is  nonlinearly 
viscoelastic. 
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The  analysis  is  carried  out  for  a  time  period  of  1200  secs.  The 
results  are  shown  in  Figs.  18-21.  At  the  free  edge  corresponding  to  x/c 
=  -1,  there  appears  to  be  the  same  kind  of  relaxation  mechanism  at  work 
since  the  stresses  at  that  location  tend  to  decrease  with  time.  Almost 
all  stresses  and  strains  exhibit  a  sharp  peak  at  x/c  =  1,  and  this  peak 
becomes  more  dominant  with  time.  The  re-entrant  corner  and  the 
associated  singularity  in  the  vicinity  of  the  right  hand  edge  is  the 
cause  for  this  unsymmetric  behavior. 

At  the  center  of  the  bond,  the  axial  and  peel  strains  do  not  change 
appreciably  with  time.  The  shear  strain  on  the  other  hand  increases  by 
24%  over  the  first  600  secs,  followed  by  a  further  increase  of  4%  over 
the  next  600  secs.  Hence  the  shear  strain  increases  asymptotically  with 
time  under  a  constant  shear  stress,  which  is  typical  for  a  viscoelastic 
material.  Numerical  integration  of  the  areas  under  the  shear  stress 
plot  at  each  time  step  yielded  the  same  value  for  the  total  shear  force 
to  within  0.5%  accuracy.  A  plot  of  the  axial  displacement  of  points 
along  the  lower  bondline  at  different  time  steps  is  shown  in  Fig.  21. 

The  displacement  of  the  upper  bondline,  shown  at  the  bottom  of  the 
figure,  does  not  change  with  time.  The  average  shear  strain  at  any 
point  along  the  bond  is  given  by  the  difference  in  the  axial 
displacement  between  upper  and  lower  bondline  divided  by  the  distance 
separating  the  two  bondlines.  The  values  of  the  average  shearing  strain 
obtained  from  the  displacement  plot  are  in  good  agreement  with  the 
results  displayed  in  the  strain  plot. 

The  stress  relaxation  mechanism  near  the  free  edge  is  not  very 
effective  at  the  right  hand  edge  due  to  the  presence  of  the 
aforementioend  singularity  which  causes  the  stresses  to  increase  with 
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time.  This  is  certainly  true  for  the  relatively  short  time  span  under 
consideration.  The  durability  of  the  bond  over  longer  time  spans  is 
likely  to  depend  upon  which  one  of  these  two  effects  dominate. 


Stress  Analysis  of  a  General  Scarf  Joint 

The  last  example  deals  with  the  stress  analysis  of  a  scarf  joint 
(see  Fig.  22a).  The  adhesive  is  FM-73  while  the  adherends  are  of 
aluminum.  The  analysis  is  restricted  to  a  constant  load  over  a  period 
of  900  sec.  under  isothermal  conditions.  The  finite  element  meshes  used 
in  these  analyses  for  a  =  0J  and  a  =  45°  are  shown  in  Figs.  22b  and  c, 
respectively. 

The  results  for  the  butt  joint  (i.e.  a  =  0°)  are  shown  in  Figs.  23 
and  24.  Stresses  and  strains  are  plotted  along  with  the  width  of  the 
joint  for  the  adhesive  and  the  adherend.  The  stresses  in  the  adhesive 
normal  to  the  bondline  are  the  same  as  those  in  the  adherend  and  their 
magnitudes  remain  constant  with  time.  Hence  axial  equilibrium  is 
satisfied  at  all  times.  On  the  other  hand,  the  normal  strain  in  the 
adhesive  is  two  orders  of  magnitude  higher  than  that  in  the  aluminum, 
and  it  shows  a  21 %  increase  in  magnitude  over  the  selected  time  span. 
Such  large  normal  strain  in  the  adhesive  gives  rise  to  fairly  large 
tensile  stress  in  the  transverse  direction,  even  though  the  transverse 
strain  is  negative.  The  shear  stresses  and  strains  are  zero  everywhere 
within  the  butt  joint,  except  near  the  free  edge. 

The  results  for  the  scarf  joint,  for  a  =  45’,  are  shown  in  Figs.  25 
and  26.  The  normal  and  transverse  stresses  follow  a  pattern  similar  to 
the  stresses  in  a  butt  joint.  The  reduction  in  the  magnitude  of  these 
stresses  is  due  to  the  45°  inclination  of  the  bond  to  the  load 
direction.  The  most  notable  difference  between  the  results  of  the  butt 
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joint  and  the  scarf  joint  is  that  the  shear  stresses  and  strains  in  the 
scarf  joint  have  large  positive  values  within  the  adhesive.  In  fact, 
the  shear  stress  is  equal  in  magnitude  to  the  normal  stress  while  the 
shear  strain  is  an  order  of  magnitude  larger  than  the  normal  strain. 

The  normal  strain  for  the  scarf  joint  increases  by  38%  over  900  seconds, 
while  the  shear  strain  increases  by  16%  over  the  same  period  of  time. 

The  large  shear  stress  present  in  the  scarf  bond  over  a  period  of 
time,  causes  the  adherends  to  slide  past  one  another.  Hence,  while  it 
may  be  advantageous  to  use  a  scarf  joint  over  a  butt  joint  due  to 
reduced  stress  levels,  long  term  loading  of  a  scarf  joint  may  lead  to 
serious  misalignment  due  to  bending  of  the  adherends. 

SUMMARY  AND  CONCLUSIONS 

A  nonlinear  viscoelastic  computational  model  is  developed, 
validated  and  applied  to  the  stress  analysis  of  adhesively  bonded 
joints.  The  nonlinear  viscoelastic  model  used  is  that  of  Schapery.  The 
finite  element  formulation  is  based  on  the  updated,  incremental 
Lagrangian  formulation.  The  program  is  validated  by  comparing  the 
present  results  with  available  analytical  and  experimental  results. 
Additional  results  for  bonded  cantilever  plate,  thick  adherend  specimen 
and  scarf  joint  are  also  presented.  In  general,  the  computer  program 
developed  herein,  called  NOVA,  is  believed  to  provide  accurate  nonlinear 
viscoelastic  analysis  capability. 

The  program  will  be  further  generalized  to  account  for  laminated 
composite  (or  anisotropic)  adherends,  moisture  effects,  and  crack 
initiation  and  growth  in  nonlinear  viscoelastic  media  in  our  future 
work . 
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TABLE  1 


Table  2 

Material  Data  for  Aluminum 
Young's  Modulus,  E^  70  x  103  MPa 

Poisson's  Ratio,  v:  0.34 

Coefficient  of  Thermal  Expansion,  a:  7.17  x  10'6  m/m/°K 
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1  Comparison  of  the  finite  element  solution  with  analytical  solution 
of  a  linear  viscoelastic  coupon  characteri zed  bv  three-narameter 


the  model  joint  problem. 
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Fig.  5  Strain  versus  time  for  FM-73  adhesive  for  step  loadinq 
at  constant  temperature. 


u.OQ  0.20  0.40  0.60  0.80  1.00  1.20  1.40  1.60 

Percent  strain 

Fig.  6  Stress  versus  strain  for  FM-73  adhesive  layer  under  linear 
temperature  and  stress. 
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Fig.  7  Stress  versus  strain  for  FM-73  adhesive  layer  under 
linear  temperature  and  constant  stress. 
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Stress,  a  (N/nmi 
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a 


x  Time  t=5  sec. 


O  Time  t=600  sec. 


qQ  =  2.5N/mm 


Distance  along  the  interface,  x  (mm) 
(a)  Axial  stress  in  aluminum  at  y=4.96mm 
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Axial  distance,  x  (mm) 

(b)  Axial  stress  in  the  adhesive  at  y=5.01  mm 

Fig.  9  Axial  stresses  in  the  adherend  and  the  adhesive  near  the 
interface  for  transverse  load  q„  =  2.5  N/mm. 
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Fig.  13  Linear  viscoelastic  plane  strain  analysis  of  a  model  joint 
(variation  of  shear  and  peel  stresses  along  theinterface) . 
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Pig.  19  Variation  of  peel  strain  and  stress  along  tne  bond 
line  for  the  thick  adherend  lap  snear  soec  *men. 
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(a)  Shear  strain  €  .  along  the  adhesive  layer 


(b)  Shear  stress  <7nt  (N/mm  )  along  the  adhesive  layer 

Fig.  26  Variation  of  shear  strain  and  stress  along  the  adhesive 
for  the  scarf  joint  (a  =45°,  aQ  =  48N/mm2). 
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